单细胞绘图数据提取个性化绘图

南山南,北秋悲

1引言

你可能有时候会不喜欢 seurat 默认绘图函数绘制的图形, 你想要好看的? 可是我怎么拿数据去画呢? 数据在哪,怎么提取呀? 我绘图技术也不咋行呀? 哎...

今天分享一些如何去提取 seurat 对象储存的数据来个性化绘图。

2数据读取

数据来自 GENES & DEVELOPMENT 单细胞结果复现 ,这里取了四个样本来测试:

library(patchwork)
library(Seurat)
library(ggplot2)
library(data.table)
library(tidyverse)
library(ggsci)
library(patchwork)
library(reshape2)

# load rds
sce <- readRDS('ythdf123_final.rds')

unique(sce$orig.ident)
# [1] ST1  ST2  ST3  ST4  ST5  ST6  ST7  ST8  SPG1 SPG2 SPG3 INT1 INT2 INT3 INT4 INT5
# [17] SER1 SER2 SER3 SER4 SER5 SER6 SER7 SER8 INT6
# 25 Levels: INT1 INT2 INT3 INT4 INT5 INT6 SER1 SER2 SER3 SER4 SER5 SER6 SER7 ... ST8

# subset data
tmp <- subset(x = sce,subset = orig.ident == c('ST1','ST2','ST3','ST4'))

unique(tmp$orig.ident)
# [1] ST1  SPG1 INT1 SER1
# 25 Levels: INT1 INT2 INT3 INT4 INT5 INT6 SER1 SER2 SER3 SER4 SER5 SER6 SER7 ... ST8

3Dimplot 函数

Dimplot 函数绘制出来的是这样的:

# dimplot
# umap
DimPlot(tmp, reduction = "umap",split.by = 'orig.ident') +
  NoLegend()

提取数据重新绘图

# get PC data
reduc <- Embeddings(tmp,reduction = 'umap') %>% data.frame()
# add group
reduc$group <- tmp$orig.ident
# add cluster
reduc$celltype <- Idents(tmp)

# check
head(reduc,3)
#                     UMAP_1    UMAP_2 group      celltype
# ST1_ATCTTGCACATC 11.644884 -1.774123   ST1 RoundSpematid
# ST1_TTTCCATAGATT 10.840847 -2.824477   ST1       Unknown
# ST1_CAGGCGTTGCTG  2.253575 12.076601   ST1  Spermatocyte

绘图:

# plot
ggplot(reduc,
       aes(x = reduc[,1],y = reduc[,2])) +
  geom_point(aes(color = celltype)) +
  theme_bw(base_size = 14) +
  scale_color_npg(name = '') +
  labs(x = 'UMAP1',y= 'UMAP2',title = '') +
  theme(strip.background = element_rect(colour = NA,fill = 'grey90'),
        aspect.ratio = 1,
        legend.position = 'top',
        plot.title = element_text(hjust = 0.5),
        panel.grid = element_blank(),
        strip.placement = 'outside',
        axis.ticks = element_blank(),
        axis.text = element_blank()) +
  facet_wrap(~group,nrow = 1)

4FeaturePlot 函数

Dimplot 函数绘制出来的是这样的:

# FeaturePlot
FeaturePlot(object = tmp,features = c("Actb","Ythdc1", "Ythdf2"),
            split.by = 'orig.ident')

提取数据重新绘图

# get gene expression
geneExp <- FetchData(object = tmp,vars = c("Actb","Ythdc1", "Ythdf2"))

# cbind
mer <- cbind(reduc,geneExp)

# merge data
megredf <- melt(mer,id.vars = colnames(reduc),
                variable.name = 'gene_name',
                value.name = 'scaledValue')

# check
head(megredf,3)
#      UMAP_1    UMAP_2 group      celltype gene_name scaledValue
# 1 11.644884 -1.774123   ST1 RoundSpematid      Actb    1.573939
# 2 10.840847 -2.824477   ST1       Unknown      Actb    1.709778
# 3  2.253575 12.076601   ST1  Spermatocyte      Actb    1.505584

绘图:

# plot
ggplot(megredf,
             aes(x = megredf[,1],y = megredf[,2])) +
  geom_point(aes(color = scaledValue)) +
  theme_bw(base_size = 14) +
  scale_color_gradient(name = '',low = 'lightgrey',high = 'red') +
  labs(x = '',y= '',title = '') +
  theme(strip.background = element_rect(colour = 'white',fill = 'grey90'),
        aspect.ratio = 1,
        legend.position = 'right',
        plot.title = element_text(hjust = 0.5),
        panel.grid = element_blank(),
        strip.placement = 'outside',
        panel.spacing = unit(0,'mm'),
        axis.ticks = element_blank(),
        axis.text = element_blank()) +
  facet_grid(gene_name~group,switch = 'y')

5结尾

只要你掌握了 ggplot 和 R 语言基础,绘图也可以很简单。

将此推文 分享到 三个微信群里 (老俊俊生信群除外), 截图发我微信, 我把代码和数据分享给你哦!

欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 (微信交流群满200人后需收取20元入群费用)。

群二维码：

老俊俊微信：

知识星球：

所以今天你学习了吗？

今天的分享就到这里了，敬请期待下一篇！

最后欢迎大家分享转发，您的点赞是对我的鼓励和肯定！

如果觉得对您帮助很大，赏杯快乐水喝喝吧！

 往期回顾 

◀...